Nintedanib induces gene expression changes in the lung of induced-rheumatoid arthritis–associated interstitial lung disease mice

Nintedanib is a multi-tyrosine kinase inhibitor widely used to treat progressive fibrosing interstitial lung diseases because it slows the reduction in forced vital capacity. However, the prognosis for patients treated with nintedanib remains poor. To improve nintedanib treatment, we examined the effects of nintedanib on gene expression in the lungs of induced-rheumatoid arthritis–associated interstitial lung disease model mice, which develop rheumatoid arthritis and subsequent pulmonary fibrosis. Using next-generation sequencing, we identified 27 upregulated and 130 downregulated genes in the lungs of these mice after treatment with nintedanib. The differentially expressed genes included mucin 5B and heat shock protein 70 family genes, which are related to interstitial lung diseases, as well as genes associated with extracellular components, particularly the myocardial architecture, suggesting unanticipated effects of nintedanib. Of the genes upregulated in the nintedanib-treated lung, expression of regulatory factor X2, which is suspected to be involved in cilia movement, and bone morphogenetic protein receptor type 2, which is involved in the pathology of pulmonary hypertension, was detected by immunohistochemistry and RNA in situ hybridization in peripheral airway epithelium and alveolar cells. Thus, the present findings indicate a set of genes whose expression alteration potentially underlies the effects of nintedanib on pulmonary fibrosis. It is expected that these findings will contribute to the development of improved nintedanib strategies for the treatment of progressive fibrosing interstitial lung diseases.


Introduction
Interstitial lung diseases are a large group of mostly progressive diseases of known or unknown causes that are characterized by chronic inflammation and fibrosis of the pulmonary interstitium. Chronic interstitial lung diseases, including idiopathic interstitial pneumonias originating from an unknown etiology and collagen vascular diseases, mostly lead to irreversible fibrotic lesions and a very poor prognosis. Although the mean survival time is 3-5 years in patients with idiopathic pulmonary fibrosis, it varies greatly among patients, making accurate assessments of prognosis difficult [1]. Although attempts to identify the causes of interstitial lung disease have been made, the pathobiological mechanisms remain elusive, leading to a range of treatment approaches. Antifibrotic drugs, nintedanib and pirfenidone, are expected to improve prognosis, and evidence has accumulated regarding their treatment efficacy in patients with interstitial lung diseases; however, in clinical practice, their efficacies vary greatly and adverse events often inhibit their treatment [2,3]. Nintedanib is a multi-tyrosine kinase inhibitor that targets platelet-derived growth factor receptor (PDGFR), vascular endothelial growth factor receptor (VEGFR), and fibroblast growth factor receptor (FGFR), and it is used clinically to slow the decrease of forced vital capacity in patients with idiopathic pulmonary fibrosis as well as in systemic scleroderma [4]. Similar effects have been reported in patients with progressive fibrosing interstitial lung disease of various origins [5]. Although the indications for nintedanib have been extended, its cost-effectiveness and frequent induction of adverse events hinder its clinical use. Thus, a better understanding of the mechanisms regulating the efficacy of and adverse events related to nintedanib treatment is needed. Anticipating that alteration of gene expression in the nintedanib-treated lung provides a clue to elucidate mechanisms regulating efficacy and adverse events of nintedanib treatment, we performed a comprehensive analysis of gene expression in the nintedanib-treated lung of a mouse model with interstitial pneumonia to identify predictors of treatment response or adverse events associated with nintedanib. Genome-wide association study is an important means of examining the relationship between gene expression and pathogenesis, and it has been used to elucidate the genetic mechanisms underlying many diseases, including idiopathic interstitial pneumonia [6]. Peyser et al. used single-cell sequencing and a bleomycin-induced pulmonary fibrosis model to characterize molecular response to fibrotic injury [7]. In that study, treatment of nintedanib reduced the bleomycin-induced cluster shift of fibroblasts to the extracellular matrix-enriched cluster, and no specific upregulation of gene expression was observed in pulmonary fibrosis lesion cells. A tissue-based sequencing approach has yet to be used to examine the gene expression alterations in nintedanib-treated lung.
Several mouse models, in which interstitial pneumonia was mostly induced by bleomycin administration via the respiratory tract, have already been used to validate the efficacy of nintedanib, including one with rheumatoid arthritis-associated lung fibrosis, in which nintedanib was found to reduce collagen levels in the lungs [8,9]. The D1CC×D1BC mouse is a doublehomozygous transgenic mouse carrying human class II major histocompatibility complex transactivator and the murine B7.1 gene under the control of the type II collagen promoter and enhancer [10,11]. Immunization of these mice with bovine type II collagen induces rheumatoid arthritis and subsequent pulmonary fibrosis in this mouse model, thereby named the induced-rheumatoid arthritis-associated interstitial lung disease (iRA-ILD) model. Histopathological and biochemical analyses have shown that these mice develop nonspecific interstitial pneumonia pattern represented by lung inflammation [12].
Here, we examined gene expression in the nintedanib-treated lung of iRA-ILD mice. We identified 27 upregulated and 130 downregulated genes in the nintedanib-treated lung and conducted functional analyses of the identified genes to elucidate their potential roles in the activity of nintedanib and the induction of adverse events. We also examined the tissue expressions of several of the identified genes. The present findings are expected to be useful for the development of predictors of treatment response or adverse events associated with nintedanib.

Mice and treatment schedule
All mouse experiments were performed according to the rules and regulations of the Fundamental Guidelines for Proper Conduct of Animal Experiments and Related Activities in Academic Research Institutions under the jurisdiction of the Ministry of Education, Culture, Sports, Science and Technology, Japan, and were approved by the Committee on the Ethics of Animal Experiments of Nagoya City University. D1CC×D1BC mice (8-12 weeks after birth and no distinction between male and female) were anesthetized with isoflurane and then immunized with 0.01 mg/ mouse of bovine type II collagen (Collagen Research Center, Japan) and an equal dose of complete Freund's adjuvant (Becton, Dickinson and Company, NJ, USA) [10,11]. The immunization time point was set at week 0, and bovine type II collagen booster injections containing incomplete Freund's adjuvant (Becton, Dickinson, and Company) were administered at weeks 3, 6, 9, and 12. The mice were then divided into two groups: nintedanib-and vehicle-treated. Nintedanib, supplied by Boehringer Ingelheim (Germany), was dissolved in 0.5% methylcellulose (Wako, Japan) and orally administered at a daily dose of 90 mg/kg from weeks 35 to 43. Blood was collected from the jugular vein for measurement of serum surfactant protein D (SP-D) levels (Yamasa, Japan). For assessment of fibrosis, the right lobe of the lung was used for the histopathological analysis, while the left lobe was used for gene expression analysis by RNA sequencing. Lungs were fixed in 4% paraformaldehyde, and 2-μm-thick paraffin sections were stained with hematoxylin and eosin and Masson's trichrome. Images showing the area of fibrosis represented in blue by Masson's trichrome staining were captured using a BZ-X analyzer (Keyence, Japan) and analyzed using ImageJ, Fiji. Data were calculated as blue ratio, with blue-stained area divided by total lung area [13]. All animal procedures, including those on D1CC×D1BC mice for in-house breeding, have been approved by the Laboratory Animal Facility of Nagoya City University.

RNA sequencing and expression profiling
For RNA extraction, samples of around 2-3 mm 3 were cut from the periphery of the frozen lungs of 3 vehicle-treated, 5 nintedanib-treated, 1 non-treated mouse, and 4 age-matched D1CC×D1BC mice, and each tissue sample was numbered. Total RNA was isolated by using an RNeasy Plus Mini Kit (Qiagen, Germany), and the integrity and purity of the total RNA were assessed by using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and a Qubit Fluorometer (Thermo Fisher Scientific). A library for RNA sequencing was prepared by using a TruSeq Stranded mRNA Preparation Kit (Illumina, CA, USA). Single reads (151 bp in length) were generated by using a MiSeq sequencer (Illumina). The reads were aligned to the mouse genome (mm10) by using TopHat [14], and computation of expression levels of annotated genes (Ensembl version 84) [15] and identification of differentially expressed genes were performed by using Cuffdiff [16]. All computational work was performed by using the DNA Data Bank of Japan supercomputer system [17].

Database analyses
Functional enrichment analysis was performed within the platform of the Database for Annotation, Visualization and Integrated Discovery (DAVID) by using categories related to biological processes, cellular components, and molecular functions [18,19]. The default value of 0.1 was used as the significance threshold in the expression analysis systematic explorer protocol. Pathway analyses were also conducted using the Kyoto Encyclopedia of Genes and Genomes and WikiPathways (http://www.wikipathways.org/index.php/WikiPathways). Each pathway was assigned a p-value, and a value of < 0.05 was considered to be statistically significant.

Quantitative reverse transcription and polymerase chain reaction (RT-qPCR)
Total RNA was extracted from the lungs as described above. cDNA was synthesized by using PrimeScript Master Mix (Perfect Real Time, Takara Bio, Japan) and amplified. Gene expression levels were compared by using an intercalator method with TB Green Premix Ex Taq II (Takara Bio) on an Applied Biosystems 7500 Fast Real-Time PCR System (Thermo Fisher Scientific). A calibration curve was constructed for comparative quantification, and glyceraldehyde-3-phosphate dehydrogenase (Gapdh) expression was used for normalization. Primer sets with the following sequences were purchased from Takara Bio (forward, reverse): regulatory factor X2 (Rfx2), GCCAGCATCACCAGCAGTACA and TCACAGTGCCTGCGATACACC; bone morphogenetic protein receptor type 2 (Bmpr2), CCATCTGAGACCTTGCTATGCTGTA and AATGTCAGC GTTCATAGTGGCATC; and Gapdh, TGTGTCCGTCGTGGATCTGA and TTGCTGTTGAAGTCG CAGGAG. RT-qPCR expression levels were compared by using Dunnett's multiple comparison test among the age-matched control, vehicle-treated, and nintedanib-treated groups. Statistical analysis was performed using the EZR software (Saitama Medical Center, Jichi Medical School, Japan), and the significance level was set at p < 0.05 [20].

Effects of nintedanib on pulmonary fibrosis and gene expression
First, we examined the effects of nintedanib on pulmonary fibrosis and gene expression in iRA-ILD mice. Histological examination revealed that nintedanib treatment alleviated destruction of the alveolar architecture ( Fig 1A: hematoxylin and eosin staining) and reduced the density of fibrotic regions (Fig 1B: Masson's trichrome staining), which is consistent with a previous report [12]. Next, using Masson staining blue ratio as an index of the level of fibrosis, we plotted the relationship between Masson stain blue ratio and serum SP-D level for 3 vehicle-treated (VA-VC) and 5 nintedanib-treated (NA-NE) iRA-ILD mice. This plot revealed a cluster of three of the nintedanib-treated mice (NA, NB, and NC) that was located far from the two vehicle-treated mice (VA and VB) ( Fig 1C). Subsequent RNA sequencing revealed 16,263 genes (Ensembl version 84) that were expressed at �1 FKPM (fragments per kilobase of transcript per million of mapped reads) in 3 or more of the 26 datasets examined (2-4 replicates derived from 5 nintedanib-and 3 vehicle-treated individuals, and 1 non-treated individual) and confirmed the presence of two major clusters: one comprising mice NA, NB, and NC, and one comprising mice VA and VB ( Fig 1D). Because mice ND, NE, and VC were not included in these two major clusters, we excluded those datasets from the following differential gene expression analysis.
A volcano plot was constructed to visualize the differential gene expression between the two clusters of mice. The results revealed significantly different gene expression for 157 genes, 27 of which were upregulated and 130 of which were downregulated (q-value < 0.05 by Cuffdiff) ( Fig 1E and Table 1). Volcano plots were also constructed to visualize the differential gene expression between the nintedanib-treated mouse with the smallest Masson stain blue ratio and each of two vehicle-treated mice (NC vs. VB, Fig 1F; and NC vs. VA, Fig 1G). Among the three combinations of mice, 3 common upregulated genes (Rfx2, Hspa1a, and Hspa1b) and 13 common downregulated genes (AC160982.1, Casq2, Fabp3-ps1, Hrc, Igkc, Igj, Myl3, Mb, Myl4, Myh6, Myom2, Tcap, and Tnnt2) were observed (Fig 2A and 2B).

Functional enrichment analysis and pathway analysis with DAVID
A functional enrichment analysis of 136 protein-coding genes from among the 157 up-or downregulated genes identified in the nintedanib-treated lungs was conducted within DAVID (Fig 3). The most significantly enriched terms were "system process" (47 genes) in the category "biological process", "contractile fiber part" (26 genes) in the category "cellular component", and "protein binding" (81 genes) in the category "molecular function". Pathway analysis using the Kyoto Encyclopedia of Genes and Genomes revealed that "cardiac muscle contraction" was the most significantly enriched pathway (10 genes) ( Table 2). In another pathway analysis tool, WikiPathways, "calcium regulation in the cardiac cell" had the highest number of six genes ( Table 2). This marked representation of cardiac functions suggested an unknown effect of nintedanib on the cardiovascular system. We therefore examined the expression of Bmpr2, which is reported to be associated with pulmonary hypertension [21] and is included among the 27 genes that were found to be significantly upregulated in the lung of nintedanib-treated iRA-ILD mice (Table 1A).

Distribution of Rfx2 and Bmpr2 expression
Next, we examined the gene and protein expression of Rfx2, one of three common upregulated genes among the three combinations of nintedanib-and vehicle-treated mice based on progression of pulmonary fibrosis (Rfx2, Hspa1a, and Hspa1b; Fig 2A), and Bmpr2. RT-qPCR and western blotting analyses showed that the expression of Rfx2 was significantly downregulated in the vehicle-treated lung compared with that in the age-matched D1CC×D1BB lung and that this downregulation tended to be ameliorated by nintedanib treatment (Fig 4A-4C). To verify the location of cells with Rfx2 expression, in situ hybridization was performed with Sftpc and Scgb1a1 transcripts, which were used as an RNA-expression marker for type II alveolar cells and for airway secreting cells, respectively. In situ hybridization revealed that Rfx2 was highly expressed in cells adjacent to Scgb1a1-positive cells localized in a basal cell-like position, and expression was observed in some type I and type II alveolar epithelial cells (Fig 4D). Subsequent immunostaining of lung tissue revealed that RFX2 was localized mainly in the bronchial epithelium, with some localization in the alveolar epithelium ( Fig 4E). RT-qPCR and western blotting analyses also showed that Bmpr2 expression tended to be downregulated in the vehicle-treated lung compared with that in the age-matched D1CC×D1BB lung, and that this downregulation also tended to be ameliorated by nintedanib treatment (Fig 5A-5C). Bmpr2 was expressed mainly in type I alveolar epithelial cells, but some expression was also observed in type II alveolar epithelial and vascular endothelial cells, as shown by in situ hybridization (Fig 5D), and bronchial epithelium and alveolar cells, as shown by immunohistochemistry ( Fig 5E).

Discussion
Using RNA sequencing, we examined the expression of genes in the lungs of nintedanibtreated iRA-ILD mice and identified 27 upregulated and 130 downregulated genes. Using functional enrichment and pathway analyses, we identified the functions and pathways enriched for these differentially expressed genes. Transcripts of two upregulated genes of interest, Rfx2 and Bmpr2, were found to be localized mainly in the bronchial epithelium and in peripheral type I alveolar cells, respectively. The RNA and protein expression levels of these genes were downregulated in iRA-ILD mouse lung compared with those in age-matched D1CC×D1BC mouse lung, and this downregulation tended to be ameliorated by nintedanib treatment.
Nintedanib is an inhibitor of multiple tyrosine kinases (i.e., PDGFR, VEGFR, and FGFR) and its effects on lung and lung fibroblasts have been validated by other groups by comparing the expressions of mesenchymal markers, cytokines, and chemokines, and by assessing changes in the expression levels of genes related to PDGFR, VEGFR, and FGFR [22]. Therefore, at the start of the present study we anticipated that we would observe changes in the expression of downstream genes associated mainly with PDGFR, VEGFR, and FGFR. However, functional enrichment and pathway analysis showed that fibrosis-related genes such as Col1a1, Fn1, and Acta2 were not enriched in the terms "signaling pathway" and "focal adhesion," which are closely related to these tyrosine kinases, indicating that nintedanib also affects the expression of genes involved in other pathways not directly associated with these tyrosine kinases. The duration of nintedanib treatment in the present study may have blurred the direct effect on the target pathway by alternative compensation or acquisition of tolerance to nintedanib-treated mice, NA, NB, and NC, and another comprising vehicle-treated mice, VA and VB. Because mice ND, NE, and VC were not included in these clusters, they were excluded from the following analysis. (E-G) Volcano plots of genes differentially expressed between nintedanib-and vehicle-treated mice. Genes with q-value < 0.05 were considered significantly up-(red dots) or downregulated (blue crosses). (E) Comparison between the two major clusters observed among the nintedanib-and vehicletreated mice. (F, G) Comparison between the mouse (NC) with the smallest Masson stain blue ratio as an index of progression of pulmonary fibrosis among the nintedanib-treated cluster and each of the vehicle-treated mice (VA and VB). The heatmap and volcano plots were drawn with the regHeatmap and plot programs in R, respectively (R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2014).
https://doi.org/10.1371/journal.pone.0270056.g001 Table 1. Genes with significantly different expression in the nintedanib-treated cluster (NA, NB, and NC) compared with that in the vehicle-treated cluster (VA and VB). inhibitors. Although comprehensive gene expression analyses have been performed in fibrotic lesions and normal regions in lungs and in lung cell lines, at the present no reports exist discussing a comprehensive analysis of gene expression in nintedanib-treated lungs [23][24][25]. Although several of the downregulated genes detected in the present study are associated with the myocardial architecture, these genes are not reported to be directly involved in the pathogenesis of pulmonary fibrosis. Although the INPULSIS trial showed a low risk of adverse cardiovascular events, it has been noted that the risk of myocardial infarction tends to increase with subsequent accumulation of nintedanib-treated patients [26,27]. On the other hand,   reports of improvements in myocardial fibrosis with nintedanib are of clinical interest [28]. The results of a recent transcriptome analysis of the mouse pulmonary myocardium with respect to atrial fibrillation have much in common with the genes downregulated in the present study [29], and the effect of nintedanib on intrapulmonary blood vessels is noteworthy. A previous study reported the effect of nintedanib on pulmonary artery smooth muscle in a rat model of pulmonary arterial hypertension [30]. Although the mechanism of how nintedanib influence genes related to cardiovascular system is not clear at the present time, the results of  , and that this downregulation tended to be attenuated by nintedanib treatment (n = 3). The vehicle-and nintedanib-treated mouse lungs were obtained from the mice used in the analysis shown in Fig 1C. (B, C) Western blot analysis confirmed the findings shown in (A): age-matched D1CC×D1BB mouse lung, n = 4; vehicle-treated mouse lung, n = 4; nintedanib-treated mouse lung, n = 4 Relative signal intensity was determined using α Tubulin as the loading control. Data are presented as mean ± standard error. � , p < 0.05 by Dunnett's test. (D) RNA in situ hybridization for Rfx2 revealed that Rfx2 was expressed strongly and weakly in cells adjacent to Scgb1a1-positive cells and in some type I and type II alveolar epithelial cells, respectively. Sftpc and Scgb1a1 transcripts were used as an RNA-expression marker for type II alveolar cells and airway secreting cells, respectively. Scale bars = 20 μm. (E) Immunohistochemical staining revealed that RFX2 (red) was localized mainly in bronchial epithelium and somewhat in alveolar epithelium. The expressions of E-cadherin (green) and SP-C (white) were used as a marker for epithelial cells, and epithelial cells and type II alveolar epithelial cells, respectively. Scale bars = 20 μm. the present functional enrichment and pathway analyses may help direct future studies. In addition, nintedanib has recently been reported to have antifibrotic activity in organs other than the lungs, including in mouse models of liver fibrosis and muscular dystrophy [31,32]; examinations for nintedanib effect in other organs are warranted.

A, Upregulated genes
In the present study, fewer upregulated genes than downregulated genes were detected, and among the 3 combinations of nintedanib-and vehicle-treated mice based on the progression of pulmonary fibrosis, only three common upregulated genes were identified (Rfx2, Hspa1a, and Hspa1b ; Fig 2A). HSPA1A and HSPA1B are members of the heat shock protein 70 family, and are closely related to one another. Several studies have reported that heat shock proteins play roles in the pathogenesis of interstitial lung disease, including in a bleomycin-induced mouse model of pulmonary fibrosis, in which induction of HSP70 expression inhibited fibrosis in the lung, and administration of anti-HSP70 antibodies reduced lung function and increased mortality [33][34][35]. In addition, extracellular matrix proteins are reported to be increased in primary lung fibroblasts lacking HSP70 [36]. In the present study, although we examined the changes in the expression levels of Hspa1a and Hspa1b in nintedanib-treated fibroblasts by means of RT-qPCR, we did not observe any changes in the expression levels of these genes at least under the present conditions.
RFX2 is a DNA-binding protein and a member of the RFX family, which is a protein family that is involved in the cell cycle, immunoregulation, and expression of motile cilia, particularly in spermatogenesis [37][38][39]. In the airway epithelium of the lung, RFX2 is suggested to coordinate with the transcription factor grainyhead-like 2 to induce differentiation of ciliated cells from progenitor basal cells during epithelial regeneration [40]. Our present immunohistochemical staining results show that RFX2 is expressed mainly in the bronchial epithelium ( Fig  5A) and our RNA in situ hybridization results show that Rfx2 is expressed in cells located adjacent to Scgb1a1-positive cells (Fig 5B), indicating that Rfx2-positive cells are localized in a basal cell-like position; therefore, nintedanib may affect airway epithelium or secretory cells via Rfx2 expression in basal cells. In addition, a recent study using single-cell RNA sequencing in pulmonary fibrosis identified pathological keratin (KRT)5 − /KRT17 + epithelial cells producing extracellular matrix and expressing the canonical basal cell transcriptional factor tumor protein p63 [41]. The Rfx2 expression-positive cells detected in the present study were located near these cells in areas of basal cells, although Krt5 showed low mRNA expression before and after nintedanib treatment. Rfx2 is anticipated to be involved in the repair of pulmonary epithelium during nintedanib treatment; however, it remains unclear how RFX2 functions in the lungs. Although we found no change in the expression of RFX family genes other than Rfx2, the regulatory genes Rfx1-3 have been analyzed by using RNA-seq and ChIP-seq in mouse ependymal cells [42]. In addition, RFX2 is reported to be involved in the regulation of the promoter of fibroblast growth factor-1B, a major transcript within the human brain and retina, by forming a complex with RFX3 in human glioblastoma cells, suggesting that RFX2 alone is insufficient for its regulation [43]. Thus, it remains unknown how nintedanib promotes Rfx2 RNA expression. To answer the question, it will be worthwhile to investigate the relationship between Rfx2 and non-coding RNAs such as microRNAs. In the present study, the expression of Mucin5b (Muc5b) was found to be significantly suppressed by nintedanib. A variant promoter of this gene, rs35705950, is reported to be a risk factor for pulmonary fibrosis; furthermore, in pulmonary fibrosis, Muc5b has been shown to be expressed not only in the conducting airways but also in epithelial cells lining honeycomb cysts, and its expression level in the mouse bronchoalveolar epithelium is associated with impaired mucociliary clearance and the extent and persistence of bleomycin-induced pulmonary fibrosis [44,45]. Taken together, the present findings indicate that nintedanib may affect the ciliary phenotypic expression and differentiation of abnormal cell populations in the peripheral airway, which may underlie its activity to slow the progression of pulmonary fibrosis.
Recently, Calabrese et al. identified 23 upregulated genes in aggregates of fibroblasts lined by epithelial cells, which they called epithelial cell/fibroblastic foci sandwich, in lungs with idiopathic pulmonary fibrosis compared to those in lungs with primary spontaneous pneumothorax by RNA sequencing analysis [46]. Intriguingly, five downregulated genes (Scgb3a1, Bpifb1, Pigr, Muc5B, and Krt5) by nintedanib treatment in the present study were included in the 14 upregulated genes with log fold-change � 2 in the 23 upregulated genes. Scgb3a1, Pigr, Bpifb1, and Muc5B are related to secretory and mucin proteins, and immune defenses [46]. Scgb3a1 is expressed in secretory cells in idiopathic pulmonary fibrosis [47]. The BPIFB1 protein was reported to be upregulated in the small airway epithelium in cystic fibrosis compared to that in the control [48]. PIGR mediates immunoglobulin A for immune exclusion of inhaled pathogens in the bronchial mucosa [49]. We previously reported Krt5 expression in hyperplastic bronchiolar cells in a mouse model of bleomycin-induced interstitial pneumonia [13]. Muc5B is reported to be expressed in the small airways and epithelial cells lining honeycomb cysts [44,45]. These findings imply that these genes could be involved in the development of pulmonary fibrosis and the effect of nintedanib in specific lesions between the peripheral airways and alveoli despite the differences in diseases, species, and sample preparation between the present study and Calabrese's study. Further investigations are needed to clarify the functions of these genes.
BMPR2 is a transmembrane and serine/threonine kinase receptor for bone morphogenetic proteins, which have functions in cell differentiation, proliferation, and bone formation [50]. Bone morphogenetic proteins are thought to affect the vascular endothelium and smooth muscle, interfering with transforming growth factor beta signaling; therefore, mutation in BMPR2 is a cause of pulmonary hypertension [21]. In addition, a previous study reported that BMPR2 expression is reduced in patients with idiopathic pulmonary fibrosis or in those with pulmonary hypertension [51]. Because we hypothesized that this blood vessel-related gene is acted upon by nintedanib via its anti-VEGFR activity, we focused on Bmpr2 among the genes upregulated by nintedanib treatment. Although Bmpr2 expression in the epithelium has been reported previously, we confirmed the presence of Bmpr2 expression in iRA-ILD mouse lung. How nintedanib upregulates Bmpr2 during the development of lung fibrosis remains to be investigated; however, we speculate that nintedanib may alleviate lung fibrosis through BMPR2 enhancement.
In the present study, Rfx2 and Bmpr2 were expressed in peripheral airway epithelium and alveolar cells, as determined by in situ hybridization and immunohistochemistry. Although it is difficult to detect cells with significant gene alteration by using current tissue-based RNA sequencing approaches, we speculate that the responsible cell population may be cells with expression of these genes. Detection of the expression of these genes is difficult in mesenchymal cells such as fibroblasts, especially in nintedanib-treated mice where the development of fibrotic lesions is attenuated by nintedanib treatment. Single-cell sequencing is expected to a useful approach for identifying the cell lineages responsible for the detected alteration of gene expression, and it is expected that it will be specific populations within the small airways or alveolar epithelium that show these changes. Single-cell sequencing has the advantage of being able to be used also for mesenchymal cell analysis, and a previous report has shown that nintedanib treatment in bleomycin-treated mice reduced a bleomycin-induced shift to the extracellular matrix-enriched cluster in fibroblasts [7], although no specific genes of activated cells were detected in pulmonary fibrosis lesions. We used fibroblasts derived from D1CC×D1BC mouse lung to evaluate the effect of nintedanib in fibroblasts; however, fibroblasts cultured after collection from lung tissue may contain several lineages, making it difficult to determine whether the established cells are the actual cells affected by nintedanib.
We acknowledge the following limitations of the present study. First, for the expression analysis, total RNA was extracted from peripheral lung tissue that included not only the alveolar epithelium, stromal tissue, pulmonary blood vessels, and peripheral airways, but also mesenchymal and hematopoietic cells such as macrophages, lymphocytes, and neutrophils. Therefore, we were unable to evaluate the differences in gene expression profiles between discrete pulmonary tissues. We did attempt to determine the specific locations of Rfx2 and Bmpr2 expression by using in situ hybridization; however, we could not quantitatively measure the changes in expression because of poor visual discrimination of the signals in specific cell types, structural differences attributed to nintedanib treatment, positioning heterogeneity in various types of hybridized cells. Localized sampling by microdissection or single-cell sequencing may overcome this limitation. In addition, although our methodology allows measurement of gene expression profiles in the periphery of nintedanib-treated lung with pulmonary fibrosis, a consistent means of distinguishing cell types will be needed to determine differential gene expression in the peripheral airway or alveoli. Second, it is difficult to conclude whether the observed up-and downregulation of genes is beneficial or detrimental to pulmonary fibrosis. Also, it remains unknown whether the observed changes in gene expression are directly involved in the antifibrotic activity of nintedanib or whether they reflect changes in response to improved pulmonary fibrosis or are the result of responses related to adverse events associated with nintedanib treatment. Additional studies based on our analysis results may lead to unexpected disadvantageous effects. Because the cellular environment of interstitial lung diseases is controlled by a complex network of interactions among mRNAs, microRNAs, long non-coding RNAs, proteins, and genomic DNA, the iRA-ILD mouse is a valuable mouse model for investigating the expression of genes in the lung environment under antifibrotic treatment.

Conclusions
Here, we used the iRA-ILD mouse model to examine the effects of nintedanib in the fibrotic lung. The present analyses revealed 157 genes that were up-or downregulated in the lung tissue of iRA-ILD mice by nintedanib treatment. Subsequent functional enrichment analysis revealed that the identified genes were associated with extracellular components, including the myocardial architecture. Of the upregulated genes, we focused on Rfx2 and Bmpr2. Rfx2 was highly expressed in cells adjacent to Scgb1a1-positive cells localized in a basal cell-like position, indicating that nintedanib treatment has some effect on the peripheral airway epithelium. Bmpr2 may be involved in the development of pulmonary fibrosis in the alveolar region. As it is difficult to collect human lung tissue for comprehensive gene expression analysis in clinical practice, this mouse model will be valuable for tracking fibrosis and treatment responsiveness over a relatively long period of medication. Our findings are expected to contribute to the development of improved strategies for the use of nintedanib for the treatment of devastating interstitial lung diseases.